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ABSTRAGT 

We study the evolution of planet-induced vortices in radially stratified disks, with initial condi¬ 
tions allowing for radial buoyancy. For this purpose we run global two-dimensional hydro dynamical 
simulations, using the PLUTO code. Planet-induced vortices are a product of the Rossby wave insta¬ 
bility (RWI) triggered in the edges of a planetary gap. In this work we assess the influence of radial 
buoyancy for the development of the vortices. We found that radial buoyancy leads to smoother 
planetary gaps, which generates weaker vortices. This effect is less pronounced for locally isothermal 
and quasi-isothermal (very small cooling rate) disks. We observed the formation of two generations of 
vortices. The first generation of vortices is formed in the outer wall of the planetary gap. The merged 
primary vortex induces accretion, depleting the mass on its orbit. This process creates a surface 
density enhancement beyond the primary vortex position. The second generation of vortices arise 
in this surface density enhancement, indicating that the bump in this region is sufficient to trigger 
the RWI. The merged secondary vortex is a promising explanation for the location of the vortex in 
the Oph IRS 48 system. Finally, we observed a nonmonotonic behavior for the vortex lifetimes as a 
function of the thermal relaxation timescale, agreeing with previous studies. The birth times of the 
secondary vortices also display a nonmonotonic behavior, which is correlated with the growth time of 
the primary vortex and its induced accretion. 

Subject headings: Accretion disks — Astrophysics: planet-disk interactions, protoplanetary disks — 
Hydrodynamics — Methods: numerical — Physics: instabilities 


1. INTRODUGTION 

High mass planets leave remarkable features in their 
parent protoplanetary disks (PPDs), namely a gap, spi¬ 
ral waves, vortices, and eccentricities. These features are 
captured i n numerical simulations of planet-disk interac¬ 
tions (e. g.. Nelson et al.|2QQQ Winters et al.||2QQ3| Klahr 




ey Dirkseh||2UU6[ |de Val-Borro et al 

2QQ71JUribe et al.|2Qll[|Lin &: PapalQizou|2Qlla|b||Atai^ 

et al. 12013 Zhu L Stone||2Q14|), and are also expected to 


et al.||2012 

IRuge et al.| 


et al.|2014t 

Juhasz et al.| 

2015 ] 


egal 

lilla 


et al. 12015 ). In this 


lution of planet-induced vortices in buoyantly unstable 
disks. 

Vortices can be formed in PPDs as a product of a 
Kelvin-Helmholtz instability, referred to as the Rossby 
wave instability (RWI) for accretion disks, and/or un¬ 
stable radial buoyancy. The RWI can be triggered when 
there is a local bu mp in the inverse potential vorticity 


profile of the disk (Lovelace et al. 1999 Li et al. 2000) 
Radial buoya ncy can be manifested as the baroclimc in¬ 


stability (BI, Klahr & Bodenheimer 2003), which needs 
a radially decreasing pressure and entropy, or in other 
words, a pressure and entropy gradients with the same 
sign. Vortices can be amplified due to the subcritical 
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baroclinic instability (SBI, Lesur & Papaloizou 2010), 
whic h is a nonlinear process. A convectiv e overstability 
(CO, |Klahr & Hubbard|j2014| |Lyra||2014| is also able to 
amplify vortices, CO is linear phase of bBI. More about 
this topic will be discussed in Section]^ Vortices are in¬ 
teresting structures to be studied, considering that they 
are important in the context of planet formation, angular 
momentum transport through the dead zone, and type I 
migration. 

In the context of planet formation, vortices are good 
candidates to trap dust grains allowing them to grow to 


planetesimal or planets size s (Barge & Sommeria 1995 
Klahr & Bodenheimer||2006|). 'This scenario is a possible 


solution for the radial drift barrier - large dust grains 
achieve high velocities toward the central star, making for 


them impossible to grow before being accreted (Whipple 
1972|). However, if the disk has a pressure bump, the 


dust grains can get trapped into this pressure maximum 


(e.g., Bracco et al.|1999| 

Varniere & Tagger|2006l Inaba & 

Barge|2UU6t|Lyra et aL|2 

009 bt Regaly et aL|2012| Meheut 

et al.||2012ap. 


gular mo mentum outwards, all o wing then matter to fall 


inwards. Shakura & Sunyaev (1973) introduced an a- 


disk model to explain this transport, where viscosity 
triggered by some kind of turbulence, is shown to be 
an efficient accretion me chanism. Usually, magne torota- 
tional instability (MRI; [Balbus fc Hawley 1991) is the 
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most invoked mechanism to explain turbulence in accre¬ 
tion disks, though in PPDs there is a region, around the 
disk’s midplane, where the level of ionization is not high 
enough for MR I to take place: the so called dead zone 


(Gammie 1996). The problem of angular momentum 
transport through t he dead zone has been investigated by 
many authors (e.g., Klahr fc Bodenheimer|2QQ3 Lesur fc 


PapaloizQu||2QlQ[ [Dzyurkevich et al.||2UlU| ' |Meh'eut et al. 


2012b|). Large-scale vortices in the dead zone of PPDs 

can hel p to transport angular momentum through that 


region. Meheut et al. (2012b) studied the angular mo¬ 
mentum urnTTaxHeToyKossby vortices. The exchange 
of angular momentum between Rossby waves in the inner 
and outer sides of a density bump, leads to a negative net 
flux, thus an outward transport of angular momentum. 

Vortices may also play a role in the context of type I 
migration. Plane t cores and l ow mass planets ex perience 
type I migration (Ward 1997 Tanaka et al.|2QQ2 ). One of 


the biggest issues about type 1 migration is the fast time 
scale in which it happens. Vortices are able to trap not 
only dust particles but also planet cores, thus t hey are 


et al. 1200311' 

Ou et al. 1120071 |Li et al.||2009| Yu et al.||20I0t 

Regaly et al 

.||20I3t jAtaiee et al.||20I4). 

4'he torm^ 

plored thorc 

rtion of planet-induced vortices is being ex- 

)ughly (e.g., Balmforth & Korycansky 20011 

de Val-Borro et al.| 20U7[ Lyra et al.| 2UU9a[ |Lin V Pa- 

paloizou||20IIb Zhu & Stone||20I4 Fu et al. 20I4| Les 


ong term 

evolution of vortices depending on the disk viscosity, 
disk temperature, and planet mass. They found criti¬ 
cal parameters for the disk viscosity (z/ = lO^^r^O^) and 
temperature {h/vp = 0.06) that lead to a long vortex 
lifetime 1 Myr). A nonmonotonic behavior with re¬ 
spect to the viscosity and temperature was found, thus 
high and low viscosities/temperatures lead to a faster 
damping of the vortices. They concluded also that disks 
with same viscosity and temperature, but more massive 
planets, in t heir case 5Mj, c an sustain vortices for a 


longer time. |Les & Lin ([2015|) studied vortex evolution 
Ls of ditterent cool 


in terms of ditterent cooling timescales. They found a 
non-monotonic dependence of the vortex lifetimes with 


the cooling timescales, which is in agreement with Fu 


et al. (2014). Moreover, they pointed out the impor¬ 
tance oil not considering locally isothermal disks, due to 
the fa ct that the RWI theory was develope d for adiabatic 
disks (Lovelace et al. |1999 Li et al.|[2QQQ ). 


In addition to the theoretical/numerical stage of this 
field, observations with high angular resolution are in¬ 
creasing. The Atacama Large Millimeter/suhmillimiter 
Array (ALMA) is now giving the capabilities to detect 
structures which may be related with unseen planets. 


ent systems: LkHa 330 (llsella et al. 2013), Oph IRS 48 

(Ivan 

der Mar el et al. 

20I3p, HD 142527 (Casassus et al. 

2013 

h'ukagawa et a 

[.I 2 OI 3 I), SAG 206462 (|Perez et al. 

2014 

), and SR 21 |P 

*erez et al. |20I4 

). An anticyclonic 


vortex could be a reasonable expt 
metrics; however, the definite explanation for these ob¬ 


servations i s still under debate (|Pinilla et al.|2015t |Flock 


The aim of this work is to study the long term evo¬ 
lution of planet-induced vortices in buoyantly unstable 


disks. The paper is laid out as follows. In Section 
we describe the planet-disk model and simulation setups. 
We describe the general evolution of our different simula¬ 
tions in Section m We discuss the role that the RWI and 
buoyancy played for vortex formation and sustenance in 
Section We study the convergence of our results with 
respect to several factors in Section We observed the 
formation of a second generation of vortices, which arise 
in a surface density enhancement that is created beyond 
the primary vortex position. The formation of the sec¬ 
ondary vortices is discussed in Section The vortex life¬ 
times and birth times with respect to different thermal 
relaxation timescales is discussed in Section Lastly, in 
Section we briefly summarize our results and state our 
conclusions. 


2. SIMULATIONS 

We study the formation and evolution of vortices in 
the outer edge of planetary gaps by solving numerically 
the following system of hydrodynamical (HD) equations 


AT 

^ + V.(Ev) = 0, 
_ + v^Vv = -—-V« 


S’ 


dp — 
_ + v^Vp. 


Sdv • V = 0, 


( 1 ) 

( 2 ) 

(3) 


where T is the gas surface density, v is the velocity, p is 
the vertically integrated pressure, is the gravitational 
potential, and Cs is the sound speed. In order to close the 
system of equations, we used an ideal equation of state 
p = c^S/y, with 7 = 1.4. 

We considered an inviscid disk, thus no prescribed vis¬ 
cosity was included. This approximation may influence 
the vortex evolution, since previous works showed that 
the vortex lifetim e is inversely proportional to the magni¬ 
tude of viscosity ( de Val-Borro et al.||2QQ7 Ataiee et al. 


2013 Fu et al. 2(JI4|). In this work, we would like to 
study the direct influence of the RWI and radial buoy¬ 
ancy for the development of the vortices. Therefore we 
chose to not include a prescribed viscosity. In our models, 
the only possible source of viscosity is the turbulence- 
triggered viscosity by the hydrodynamical instabilities. 
Lastly, we assumed that the barycenter of the system is 
located at the star’s center. This simplification is plau¬ 
sible, because the planet masses considered are not very 
large (I and 3Mj) neither the vortices accumulate much 
massthus the deviation of the barycenter with respect 
to the star’s center should be small. Nonetheless, this ap¬ 
proximation may slightly influence planet-induced vortex 
formation, since it eliminates the Lagrange point L3, in 
the corotation region, which could change the gap struc¬ 
ture. 

We used the pla net-disk module fo r the PLUTO code 
that is presented in Uribe et al. (2011). The gravitational 
potential includes contributions froni the planet and the 
star, and it is given by 




GM, 


P 


GM, 


(4) 


y(|r-Rp|2+e2) M’ 

^ We obtained vortices masses up to a few 10“^ Mq, integrating 
the surface density with respect to the area element inside the 
vortex region. 
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where G is the gravitational constant, Mp is the planet 
mass, Rp is the planet location, is the stellar mass 
and e is a softening parameter. It is needed to soften the 
gravitational potential of the planet in order to avoid nu¬ 
merical divergence close to the planet’s location. More¬ 
over, this softening can account for 3D effects of vertical 
stratification. We considered this parameter as being a 
fraction of the Hill radius e = kRn, with k = 0.6 


8ind Rh = Rp{Mp/. The recommended soft¬ 
ening factor for t he planet gravita tional potential is of 
e = 0.6i^ — d.lH ( |Kley et aL||2Q12 ), where H is the disk 
scale height. These values can recover 3D effects of ver¬ 
tical stratification. The Hill radius and the disk scale 
height at the planet position are similar in our simula¬ 
tions, thus we chose e = d.dRn- 

The stationary solution of a sub-Keplerian disk was 
taken as initial conditions, which in polar coordinates is 
given by 


E = Eo 


Cs = Co 



Vr = 0 , 


(5) 

( 6 ) 
(7) 





r dp 


( 8 ) 


where Eq is the initial surface density at ro = 1 AU, 
/3 e = 1.5 is the slope for the power law distribution of 
surface density, = 0.5 is the slope for the power law 
distribution of sound speed, vk is the Keplerian velocity, 
and h = CsjvK — Hir = 0.05 is the initial aspect ratio 
and fix the intial sound speed cq at tq = 1 AU, since vk 
at ro = 1 AU is set as one. 

The planet is set up as a point mass in a given position 
Rp and with a given mass Mp . In order to avoid an initial 
big disturbance to the disk, we added the planet slowly 
along its first Keplerian orbit, according to the following 


M'p = Mp 



where t is the global time and P = 27r{GM^lRp)~^^‘^ is 
one planetary orbit. Thus, while t < T the planet mass 
slowly increases toward Mp. 

The planet initial velocity is assumed to be the Keple¬ 
rian velocity = yjGM^JWp and the initial accelera¬ 
tion coming from the gravitational interaction with the 
star and the disk is 


ap — 


GM.Rp 

|Rp P 



GI](r - Rp) 
y(|r-Rp|2 +e2)3 


dA, (10) 


where dA is the area element and ^ is a factor that soften 
the contribution of the disk gravity in the Hill sphere and 
is given by 


^ = 1 — exp 


|r-RpP l 

(0.6i?ff)2 J- 


( 11 ) 


Lastly, the planet position, velocity, and acceleration 
changed according to the dynamical interaction with the 
star-disk system. Its acceleration changes with time fol¬ 
lowing Equation \To \ and its position and velocity are then 
updated using a leapfrog integrator. 


2.1. Thermal Relaxation 

In order to account for radiative effects, we applied 
cooling to the system. We modeled this cooling via ther¬ 
mal relaxation, following the approach below 

dT _ (T - T^) 
dt r(r) 

where T is the temperature, is the initial temperature 
(equilibrium temperature as result of irradiation), and 
r(r) is the relaxation timescale, which depends on radius 
(r(r) = 27rr/U(r)). This approach tends to reestablish 
the equilibrium temperature profile, after the planetary 
gap is opened and the system reaches a steady state. 

Instead of adding cooling as a source term in the energy 
equation, we updated the temperature at each time step 
according to equation Numerically it corresponds to 

2^new _ 2^old _ _ T^) (f8) 

where is the relaxed temperature, would be 

the temperature we get from the solution of the energy 
equation, and At is the time step. We solve Equation 
which describes conservation of energy, considering 
pressure as an independent variable. Hence, we had to 
convert equationfrom temperature to pressure depen¬ 
dent, for which we used the relation T (x p/E, leading 
to 

where p^^^ is the new pressure from the relaxed temper¬ 
ature, is the pressure we get from the solution of 
equation p^ is the initial pressure, E^®^ is the den¬ 
sity we get from the solution of equation |h and E^ is 
the initial density. Einally, we cooled the disk through 
equationEor the locally isothermal setups, instead of 
using Equation \T^ to cool the disk, we setup the sound 
speed to its initial profile at every time step, in order to 
guarantee locally isothermality. 

2.2. Numerieal Setup 

The simulations were carried ou t using the finite vol - 
ume Godunov-type code PLUTO ( [Mignone et al.|2QQ7 ). 
Spacial integration and time evolution were performed 
using the piecewise parabolic method (PPM) and second 
order Runge-Kutta schemes, respectively. The Harten- 
Lax-van Leer-Contact (HLLC) Riemann solver was used 
to compute the numerical fluxes and the Strang oper¬ 
ator splitting method to solve the equations in multi¬ 
dimensions. 

The HD equations were solved in a two-dimensional do¬ 
main considering polar coordinates. A logarithmic grid 
was used for the radius and a uniform one for azimuth. 
The system was integrated from 0.25 AU to 4.0 AU in 
radius and from 0 to 27r in azimuth. Temporal evolu¬ 
tion was taken up to 5000 orbits. Reflective boundary 
conditions were used in the radial direction and peri¬ 
odic conditions in the azimuthal direction. Distances 
are given in units of 1 AU; surface densities in units of 
Eo = 10“"^ Mq/AU^, which corresponds to a disk mass 
of O.OO 2 M 0 inside the domain considered, therefore it is 
plausible to neglect disk self-gravity, since the Toomre 
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parameter is Q >> 1 everywhere in the disk; and ve¬ 
locities in units of Keplerian speed at 1 AU. Table 
summarizes the simulations parameters. 


Table 1 

Simulations parameters 


Label 

Mp^ 

(Mj) 

Rp^ 

(AU) 

{2tt/Q,o) 

{Nr,N4,)‘^ 

TROOl 

1.0 

1.0 

0.01 

(512, 1024) 

TROl 

1.0 

1.0 

0.1 

(512, 1024) 

TRl 

1.0 

1.0 

1.0 

(512, 1024) 

TR2 

1.0 

1.0 

2.0 

(512, 1024) 

TR5 

1.0 

1.0 

5.0 

(512, 1024) 

TRIO 

1.0 

1.0 

10.0 

(512, 1024) 

ISOIMJ 

1.0 

1.0 

0.0 

(512, 1024) 

IS03MJ 

3.0 

1.0 

0.0 

(512, 1024) 


^ Planet mass in terms of Jupiter mass (consider¬ 
ing M* = Mq). 

^ Planet location in AU. 

Thermal relaxation timescale in orbital units. 

^ Numerical resolution in the radial (Nr) and az¬ 
imuthal (A^^) directions. 


3. GENERAL EVOLUTION 

In this section we describe the system evolution for 
our simulations. Eirst, we present the results for the 
simulations with a IMj planet mass, varying the thermal 
relaxation timescales. Second, we present the results for 
the isothermal simulations, considering IMj and 3Mj 
planet masses. 


3.1. Non-isothermal Cases 

The first set of results presented are the cases with a 
IMj planet and different thermal relaxation timescales. 
The different values of considered and their labels 
can be seen in Table All the simulations presented a 
similar behavior, whidi we describe hereafter. 

The formation of spiral waves takes place during the 
first planetary orbit. Additionally, during the first tens 
of orbits the planet carves out a very noticeable gap and 
small vortices are formed in the outer edge of this gap. In 
the first hundreds of planetary orbits these small vortices 
merge into a bigger one. Some mass remains in the gap 
region even after a few thousands of orbits, which can be 
related to the inviscid and/or non-barycentric approxi¬ 
mations. The inviscid approximation may influence the 
efficiency of mass transport. Nonetheless, neglecting the 
indirect potential exerted on the disk due to the barycen- 
ter shift also seems to retain mass in the gap region, even 
for non-inviscid disks (see Eigure 2 in Zhu & Stone|2014). 
In the first thousands of planetary orbits a surface den- 
sity enhancement appears beyond the vortex position. 
Accumulation of mass persists and a second vortex is 
formed in this region. The primary vortex is damped 
in different timescales for the different Ur’s, nonetheless 
some material remains in a ring-shape form in between 
the planetary gap and the secondary vortex. This mate¬ 
rial is also dispersed out in different timescales. 

Simulations TROOl and TROl present a secondary vor¬ 
tex very similar to the primary one. Eigure shows the 
system evolution for TROl. Simulations TRl, TR2, TR5, 


® Hereafter, we refer to r(r) as Q.t. 


and TRIO present also a secondary vortex; however, the 
new vortex is highly spread in the azimuthal direction. 
EigureJ^ shows the potential vorticity at 5000 orbits for 
the different Ur’s. For Ur = 0.01, the secondary vor¬ 
tex survives until the end of the simulation; however, the 
vortex gets spread along the azimuthal direction. For 
Ur = 0.1, the secondary vortex survives and does not 
get spread in the azimuthal direction. For Ur > 1.0, 
the secondary vortex is mostly damped by the end of the 
simulation interval. 


3.2. Isothermal Cases 

Here we present the results for the isothermal setup 
and planet masses of IMj and 3Mj. The simulations 
labels are presented on Table 

The isothermal configuration shows a considerably 
similar behavior as the models with thermal relaxation. 
For the ISO IMJ simulation, the sequence of events is the 
same. We first observe the formation of spiral waves, fol¬ 
lowed by planet gap opening, and production of small 
vortices at the outer edge of this gap. The small vortices 
gather together and merge into a bigger one. A surface 
density enhancement appears beyond the primary vor¬ 
tex position. Material is accumulated at this location 
and a second vortex arises, this structure gets spread 
in the azimuthal direction with time. The primary vor¬ 
tex gets damped and the material in between the plan¬ 
etary gap and the secondary vortex disperses out. The 
timescales for the events are similar to the ones for the 
non-isothermal cases. 

The IS03MJ simulation presents a similar sequence of 
events, with the difference that two vortices, that do not 
merge with time, are created in the outer edge of the 
planetary gap. The evolutionary timescales for which 
different structures form are also different. The surface 
density enhancement appears in hundreds of planetary 
orbits, instead of thousands. The damping of the pri¬ 
mary vortices is also faster. Pile-up of material at the 
surface density enhancement also happens. Nonetheless, 
it takes thousands of planetary orbits for a secondary 
vortex to arise. After a few thousands of planetary or¬ 
bits the material between the planetary gap and the sec¬ 
ondary vortex is totally dispersed out, and a much wider 
gap is settled. Figure shows the system evolution for 
simulation IS03MJ. 

It was not possible to consider a higher planet mass 
(e.g., lOMj) under the setup assumed, because the gap 
created is much wider, which makes the disk size consid¬ 
ered too small. To solve the same structures in a bigger 
disk, we would need to use more grid cells. 


4. VORTEX FORMATION AND EVOLUTION 
In a protoplanetary disk, we know that vortex forma¬ 


tion can happen as a product of the RWI (Lovel ace et al 
1999| Li et al.l 20001) and/or radial buoyancy ( Klahr V 


Bodenheimerl. 

Hubbard 2Q14| ) 


l[2QQQt _____ 

2UU3[ [Lesur fc Papaloizon||2QlQ| Klahr V 
iT In this work, we considered initially 


buoyantly unstable disks. Nonetheless, we know that the 
presence of a planetary gap naturally triggers the RWI, 
due to the sharp surface density gradient that is created 
in the gap edges. Hereafter, we discuss the role that the 
RWI and radial buoyancy played for the formation and 
evolution of planet-induced vortices. We would like to 
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Figure 1. Evolution of the surface density perturbation (top panel) and the potential vorticity with the Keplerian profile subtracted 
(bottom panel). The color bar for the potential vorticity plots was truncated from —0.5 to 0.5, in order to provide a higher contrast. The 
results show simulation TROl. 



Figure 2. Final potential vorticity with the Keplerian profile subtracted for the different f^r’s considered. 


remember that we are using an inviscid approximation, 
thus any viscosity in the system is turbulence-triggered 
viscosity by the hydro dynamical instabilities. We chose 
to consider an inviscid approximation to assess the direct 
influence of radial buoyancy and the RWI for the vortices 
evolution. 


4.1. Rossby Wave Instability 

The RWI is a pressure driven instability for rotating 
systems, which is composed of non-axisymmetric modes. 
The RWI is triggered when there is a local maximum in 
the radial profile of the function ( [Lovelace et al.|p!9^ ) 


C{r) = 


(15) 
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Figure 3. Evolution of the surface density perturbation (top panel) and the potential vorticity with the Keplerian profile subtracted 
(bottom panel). The color bar for the potential vorticity plots was truncated from —0.5 to 0.5, in order to provide a higher contrast. The 
results shows simulation IS03MJ. 


where = (V x fT) • z/H is the potential vorticit 3 |^and 
S = p/TH is an equivalent to the entropy. Physically, this 
condition can be achieve d at the edge of planetary gaps 
(de Val-Borro et al.|2007) and the edge of dead zones d ue 
to sharp viscosity transitions (Lyra & Mac Low|2012). 

The formation of vortices as a product of the KWI 
has been studied by many authors. The gro wth rate 
of the inst ability for 2D dis ks was studied by |Li et al. 
(2000) and Umurhan (2010), usin g different approxima¬ 


tion s, and in 3D stratified disks by 


and Lin 


was exp 


(2012b ). The nonlinear p 


Meheut et al. (2012b) 
lase of the instability 


ored by Meheut et al. (2013). Despite t he theory 
for th e RWI was developed lor adiabatic disks (Li et al. 
2000), most of the works used locally isothernial disks 
to study the formation and evolution of planet-induced 
vortices (Balmforth fc Korycansky 200ll de Val-Borro 


_ Jxoryi 

et al.||200V |Lyra et^ir|2009al |Lin fc Papaloizou||201lR 


Zhu et al. I 20141 |Fu et al.||2014p. Just recently, |Les & 

Lin| ( |2015| ) made a breakthrough and added an artiti- 


cial source of cooling and heating to explore the non- 
isothermal behavior. In this context, our work is a second 
step to the process of understanding the non-isothermal 
scenario. 

It is not in the scope of this work to make an extensive 
study of the growth and decay of the RW I , sinc e this 
matter was already addressed by Les & Lin (2015) for a 
physical scenario very similar to ours. In order to have 
just a qualitative insight, we analyzed the spacetime evo¬ 
lution of the potential vorticity (C) averaged in azimuth. 
Figured shows the result for simulation TROl. 

The blue region in Figure represents a minimum in 
(■. A minimum in C is equivalent to a maximum in C 
(Equation [T^, sufficient condition to trigger the RWI. 


^ Hereafter, called as C instead oi T ^. 


a 



Figure 4. Spacetime evolution of the potential vorticity (C) av¬ 
eraged in azimuth for simulation TROl. The Keplerian profile was 
subtracted. The color bar was truncated from —1.0 to 1.0 in order 
to obtain a higher contrast. 


This minimum is achieved in both, the outer edge of the 
planetary gap and at the surface density enhancement 
outwards the primary vortex position. Therefore, the 
RWI has been triggered in both regions. The presence 
of the minimum is maintained along the whole simula¬ 
tion interval. In the planetary gap edge, its value slowly 
increases with time, which explains the vortex decay. In 
the surface density enhancement outwards the primary 
vortex position, its value is kept slightly constant with 
time, which explains the survival of the secondary vortex 
until the end of our simulation. The spacetime evolu¬ 
tion of C is similar for all the cases. A local minimum 
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is observed in the regions of the primary and secondary 
vortices. The main differences are the size of the blue 
regions (local minimums), the time the local minimum 
related the primary vortex starts to decay, and the time 
the local minimum related to the secondary vortex ap¬ 
pears. 


4.2. Radial Buoyancy 

The necessary ingredients for a CO and SBI are: (i) ra¬ 
dial pressure and entropy gradients possessing the same 
sign (radial buoyancy) and (ii) thermal relaxation, with 
maximum amplification for Dr 1.0. The f ormation 
of vortices due to th e BI was first obser ved by|Klahr fc: 
Bodenhei merl (|2003|) . Further studies by |Petersen et al. 
(|2UU7a|d ) showed tne importance of thermal relaxation 


for baroclinic vortex amplification. They found that 
thermal relaxation or diffusion, besides entropy gradient, 
are required to k eep the instability in action. [Lesur fc 
Papaloizouj (|2010|) studied baroclinic vortex amplihcation 
through the growth of existing vortical perturbations. In 
order to not cause confusion between the generation of 
vortices by the classical BI and amplification of the vor¬ 
tices in a radial buoyant fluid, they coined this process a 
SBI. A parametric study covering the important ranges 
of entropy gradients, thermal diffusion timescales, and 
thermal relaxation timescales for PPDs was carried out 


by Raettig et al. (2013). They showed the importance of 
baroclinic ettects even for small entropy gradien ts, which 


is the case in PPDs. Klahr & Hubbard (2014) found a 
linear amplification oil epicyclic oscillations m radially 
stratified and vertically unstratified disks, which they 
called convective over stability. This phenomenon can be 
regarded as the linear phase of the SBI. Yet not much ef¬ 
forts were made to study how vortex formation and am¬ 
plification proceeds in a buo yantly unstable d isk with a 
high mass planet embedded. Les & Lin (|2015 ) discussed 
briefly whether an axisymmetric instability was at play 
in their simulations of planet induced vortices; however, 
they concluded that only the RWI was in action. 

We can quantify the radial stability in a disk with re¬ 
gards to convection thro ugh the Brunt-Vaisa la frequency 
(V), which is given by ( Raettig et al.||2Q13 ) 


= - 




(16) 


where /3p is the pressure gradient, Ps is the entropy gra¬ 
dient, and D is the angular velocity. Positive values of 
N‘^ indicate stability. The entropy gradient for a 2D ver¬ 
tically integrated disk is given by 


/^S' = /^T + (1 - 7)^E, 


(17) 


where Pt is the temperature gradient and is the sur¬ 
face density gradient. We made the choice for the initial 
surface density and sound speed gradients in a way that 
it gives an initial negative value for N‘^ equals to —0.0018, 
thus favoring instability. 

Figure]^ shows the spacetime evolution of N‘^/Q? av¬ 
eraged in azimuth for simulation TROI. We plot 
instead of Y^, to eliminate the dependence with the an¬ 
gular velocity. Since Q? is always positive, Y^/D^ > 0 
still indicates stability. We can see that in the outer disk 
Y^ is kept negative and roughly equals to its initial value, 
with exception for the outer boundary. The outer radial 


extent with negative Y^ becomes narrower throughout 
the simulation interval, thus the evolution of the system 
tends to stabilize the disk with respect to buoyancy. In 
the gap region and outer gap wall, Y^ becomes positive 
after a few tens of planetary orbits; however, there is 
a strip around the primary vortex position with smaller 
values of Y^. The strip’s center possesses negative Y^ 
in the first tens of planetary orbits, but Y^ turns posi¬ 
tive later on. From ^ 600 orbits, Y^ becomes negative 
again in the central position of the vortex. The region 
with negative Y^ is not as large as the vortex size, thus 
buoyancy is not playing a major role for the evolution of 
the primary vortex. Y^ is positive in the region where 
the second generation of vortex appears, until the time 
that the secondary vortex arises and a strip with negative 
Y^ appears around the secondary vortex position. Once 
more, the strip width is smaller than the vortex size, in¬ 
dicating that buoyancy may not be playing a major role 
in the location of the vortex. 

The behavior of Y^ for the other cases is similar to 
the one of TROI, with exception for TROOI and the 
isothermal cases. For them, Y^ becomes positive in the 
first tens of planetary orbits and remains positive along 
the whole simulation interval. This indicates that buoy¬ 
ancy does not play any role for the quasi-isothe rmal and 
isothermal case s. Reinforcing the findings of Petersen] 
et al. ( |2007a|b| |, in the lack of thermal relaxation or dif¬ 


fusion, buoyancy is not sustained. For Dr > I.O, the 
strip around the position of the primary vortex has neg¬ 
ative Y^ during a larger fraction of the primary vor¬ 
tex lifetime. The strip around the secondary vortex is 
wider, indicating that buoyancy had more importance 
for the secondary vortex evolution than in the case that 
Dr < I.O. 
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Figure 5. Spacetime evolution of N‘^/Q? averaged in azimuth for 
simulation TROI. The color bar was truncated from /Q?‘)ini 

to in order to provide a higher contrast. 


To check the impact that buoyancy has in the results, 
we used a model where the initial Y^ is positive. We run 
the new simulation with the same physical and numeri¬ 
cal setup as for the TROI case, but changing the surface 
density gradient from = 1.5 to = 3.0. The gen¬ 
eral evolution of the system was very similar to the case 
where Y^ is initially negative. Two major differences 
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were noticed. The first is regarding the maximum ampli¬ 
tudes that the primary and secondary vortices achieve, 
which is higher for the case where N‘^ is initially posi¬ 
tive. The second is regarding the second generation of 
vortices. For the > 0 case, two vortices arise in the 
surface density enhancement region and take more time 
to merge 1000 orbits against ^ 500 orbits for the ref¬ 
erence case). The secondary merged vortex has also an 
aspect ratio much smaller than the secondary vortex for 
the Nl, < 0 case. Figure [^presents a comparison of the 
surface density profiles for the positive and negative 
cases, for two different points in time (t = 250 orbits and 
t = 2500 orbits). 


t = 250 orbits 



r(au) 


Figure 6. Surface density profiles averaged in azimuth. The red 
dashed line shows N‘^ initially negative, whereas the slate blue 
dotted-dashed line initially positive. 

The planetary gap structure is very similar for the dif¬ 
ferent surface density slopes, the same width and lower 
level for the depth are observed, as well as the same 
location for the maximum and minimum surface den¬ 
sity perturbations. The standing difference is the sharp¬ 
ness of the surface density gradient in the planetary gap 
edge and at the surface density enhancement beyond 
the primary vortex position. A larger gradient produces 
stronger vortices, therefore the vortices for > 0 are 
stronger. The times chosen to compare the cases rep¬ 
resent a moment the vortices are totally merged. The 
spacetime evolution of N‘^/Q? for = 3.0 (Figure]^ 
shows that N‘^ becomes negative in the first tens of orbits 
in the region where the primary vortices arise; however, 
it becomes positive again and local buoyancy disappears 
for hundreds of orbits, N‘^ becomes negative again from 
^ 700 orbits. The strip with negative N‘^ around the 
vortex position is again not as wide as the vortex size, 
indicating that buoyancy is not playing a major role for 
the evolution of the primary vortex. Nonetheless, this 
shows that initially buoyantly stable disks can undergo 
an inversion of sign for the entropy gradient, therefore 


turning on instability. In the secondary vortex region, 
N‘^ never turns to be negative. The RWI is the only 
responsible for the formation and sustenance of the sec¬ 
ondary vortex. 



r(au) 

Figure 7. Spacetime evolution of N‘^/Q? averaged in azimuth for 
the simulation with > 0. The color bar was truncated from 

— {N‘^/Q?)ini to /Q?)ini^ in order to provide a higher contrast. 

This result demonstrates that when we have the RWI 
and buoyancy acting at the same time, weaker vortices 
are produced. Therefore, buoyancy opposes vortex am¬ 
plification and survival, in this scenario. Taking into 
account that these processes provide viscosity to the sys¬ 
tem, once we have both in action more viscosity is pro¬ 
duced. More viscous disks carve out smoother gaps, lead¬ 
ing to the weaker vortices. It should also be noticed that 
for Ur > 1.0, the secondary vortices get damped dur¬ 
ing the simulation interval, those are also the cases for 
which buoyancy plays some role for the secondary vortex 
evolution. 

5. CODE CONTROL 

In this section, we explore the numerical factors that 
could influence the physical validity of our simulations. 
The tests were done using the TROl case as reference. 
The physical conditions and numerical setup were ex¬ 
actly the same as for TROl, varying only what we follow¬ 
ing mention. We checked the convergence of the results 
considering a higher numerical resolution. We analyzed 
whether our reflecting boundary conditions for the radial 
direction may have reflected waves, influencing the evolu¬ 
tion of the system. A different way to prevent boundary 
effects is to push the outer disk to a larger radius, hence 
we also used this approach to check whether the disk size 
influenced the results. Lastly, we added the planet along 
a larger number of planetary orbits to analyse whether 
the initial planet disturbance could have been too large, 
generating fake effects. 

For the numerical resolution test, we doubled the reso¬ 
lution from N(f)) = (512, 1024) to (A^^, A^^) = (1024, 
2048). The temporal evolution was taken up to 1000 
planetary orbits. The full temporal evolution was not 
checked, because the doubled resolution is numerically 
highly expensive. For the boundary conditions test, we 
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first changed the inner and outer radial boundary con¬ 
ditions from reflective to non-reflective, second we con¬ 
sidered a larger disk extending from 0.25 AU to 8 AU. 
Finally, aiming to check the effect of the initial planet dis¬ 
turbance to the disk, we made two tests: slowly adding 
the planet along its first 10 and 100 orbits, following 
Equation in the reference case the planet was slowly 
added along its first orbit. The temporal evolution for 
the last four tests was taken up to 5000 planetary orbits. 

We compare our test cases with the reference case 
using their surface density profiles at the latest snap¬ 
shot. Figure [^presents these profiles for t = 1000 orbits 
(numerical resolution comparison) and t = 5000 orbits 
(other comparisons). We observed a good agreement for 
the surface density profiles, indicating that the simula¬ 
tion results were not much influenced by these factors. 
Nonetheless, for the resolution test, the outer gradient 
that leads to the second generation of vortices is slightly 
smoother than for the standard case. 

For the boundary condition test, in the outer disk 
(r > 3 AU) the material is emptied out, due to the 
boundary choice; however, the main results agree with 
the standard case. The major difference is regarding 
the secondary vortex, for the standard case a secondary 
merged vortex is created, in this case two secondary vor¬ 
tices are created and they do not merge until the end of 
the simulation. The two secondary vortices, for the non¬ 
reflecting radial boundaries, are right opposite to each 
other in the azimuthal direction and they have about the 
same strength. We speculate that the flow in the coro¬ 
tation region of these vortices is not being able to push 
them together, exactly because the vortices have same 
strength and are located right opposite to each other, 
leading to a stable configuration. For the larger disk size, 
the standing difference is regarding the local minimum 
between the surface density bumps, where the primary 
and secondary vortices sit. For the standard case, the 
local minimum is still present by the end of the simu¬ 
lation. For the larger disk size, this local minimum has 
disappeared by the end of the simulation, meaning that 
the matter is dissipating slightly faster. 

For the planet being added along the first 10 orbits, 
the primary small vortices merge in a faster timescale, 
and the vortices amplitudes take a longer time to grow 
than for the standard case. For the planet being added 
along the first 100 orbits, the vortices amplitudes take 
also a longer time to grow than for the standard case. 
Higher pla net masses can excite stronger vortices ( |Fu| 
et al. 120141 ), thus when the planet perturbation is addeo 


slower, the vortices will also take a longer time to grow. 


6. SECOND GENERATION OF VORTICES 

The most interesting result of our simulations is the 
second generation of vortices. A surface density enhance¬ 
ment was observed beyond the primary vortex position 
for all the cases. This bump is strong enough to trigger 
the RWI outside the primary vortex radius and to form 
a second generation of vortices. As it was discussed be¬ 
fore, we observe a minimum for C in the region of the 
secondary vortex for all the cases. Also, there are strips 
of negative N‘^ in the region of the secondary vortex, with 
exception for the TROOl and the isothermal cases. Once 
more, the action of the RWI together with buoyancy con¬ 
trols the vortex evolution. For the TROOl and isother¬ 



0 12 3 4 

r (au) 

Figure 8. Surface density profiles (averaged in azimuth) at 
1000 orbits (top panel) and 5000 orbits (bottom panel). The top 
panel presents the comparison for the numerical resolution test. 
The red solid line shows the reference case, whereas the slate 
blue dotted-dashed line the numerical resolution test. The bot¬ 
tom panel presents the comparison for the other cases. The red 
solid line shows the reference case, the goldenrod dotted line shows 
the boundary conditions test, the green dashed line shows the disk 
size test, the slate blue dotted-dashed line shows the planet distur¬ 
bance test (the planet was added during its first 10 orbits), and the 
violet triple-dotted-dashed line shows the planet disturbance test 
(the planet was added during its first 100 orbits). 

mal cases, the RWI is solely responsible for the secondary 
vortex. We already established that our results were not 
affected by the choice of boundary conditions, resolution, 
or the planet perturbation being too sharp. Therefore, it 
confirms the physical origin of the second generation of 
vortices. No further density enhancement (strong enough 
to keep triggering the RWI) beyond the secondary vortex 
position was observed, even for the test simulation with 
a larger disk, thus we do not expect a third generation 
of vortices. 

6.1. The Origin of the Seeondary Vortex 

Accumulation of mass is observed at the inner bound¬ 
ary for all the cases. Our understanding is that the 
primary vortices produce an effective a-viscosity that 
is large enough to induce accretion in the disk in the 
timescales we simulate. Therefore mass is flowing from 
the region of the primary vortex to the inner disk. The 
depletion of mass in the orbit of the primary vortex po¬ 
sition looks like a gap carved out by the primary vortex. 
In fact, the outer wall of the planetary gap is moving 
outwards due to this depletion of mass. For instance, for 
the IS03MJ case, by the end of our simulation all the 
mass in the orbit of the primary vortex was already de¬ 
pleted and the planetary gap became wider in its outer 
side. Figure [^presents the inner disk mass as a function 
of time for simulation TROl, where the increase of the 
mass is demonstrated. 

The faster mass increase happens in the first tens of 
orbits, when the planet still carving out its gap. During 
this period, the mass increase is mostly due to the plan- 
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Figure 9. Inner disk mass as a function of time for simulation 
TROl. We consider 0.25 AU < r < 0.75 AU as the inner disk. 


etary gap opening process. The inner disk mass is being 
pushed from r = 1 AU to the inner parts of the disk. 
For t > 100 orbits, the mass increase is most likely due 
to accretion induced by the primary vortex. Our simu¬ 
lations assumed an inviscid disk, therefore any viscosity 
is produced by the hydro dynamical instabilities, which 
in this case could be named a vortex-induced viscosity. 
We can obtain an estimate of the ^-parameter, using the 
r — (f) compo nent of the Reyno lds stress and the local 
sound speed (|Flock et al. 2011|). For simulation TROl, 
the Qf-parameter averaged m space and time (until the 
time the secondary vortex appears) was a ^ 3 x 10“^, 
a value in the range of wha t is typically obtained by 
the MRI, a = 10“^ — 10“^ (Dzyurkevich et al. 2010). 
The other simulations presented values in the range from 
a = 10“^ — 10“^, in agreement with MRI. The lowest 
value obtained was a 6 x 10“^, for TRl, and the high¬ 
est value was a 1 x 10“^, for TR5. Large-scale vor¬ 
tices are able to transport angular momentum outwards, 
because a negative angular momentum flux is obtained 
from the balance between the angular momentum car¬ 
ried by Rossby waves in the inner and out er sides of a 
surface density bump (Meheut et al.||2012b). 

In this work we assumed that the barycenter of the 
system is located at the star’s center. This approxima¬ 
tion could influence the gap structure, since the Lagrange 
point L3, in the corotation region of the planet, is re¬ 
moved. Changes in the gap structure may affect the pri¬ 
mary vortex generation, subsequently possibly impact¬ 
ing the second generation of vortices. This assumption 
summed to the inviscid disk approach are two factors 
that could influence the formation of a second genera¬ 
tion of vortices. We suggest that further studies should 
check these factors. 


6.2. Pressure Bumps 

Vortices are able to trap dust particles, because they 
are local pressure bumps. The particles are attracted to 
the highest pressure region, thus to the vortex center. 
The secondary vortex becomes extensively spread in the 
azimuthal direction for Ur >1.0 and Mp = 3Mj. In or¬ 
der to make sure that these nonaxisymmetric structures 
can still trap dust particles, we checked their pressure 


profiles. Figure pR| shows a radial cut of the pressure for 
equals to the vortex center and an azimuthal cut of the 
pressure for r equals to the vortex center. We show the 
cases of Ur = 10.0 and Mp = 3Mj, since these are the 
simulations that present the most spread vortices in the 
azimuthal direction. We can observe a pressure bump in 
both radial and azimuthal directions; however, in the az- 
imnthal direct i on th e bump of IS03MJ is very smooth. 


Birnstiel et al. (2013) showed that a very smooth pressure 
bump in the azimutnal direction is still sufficient to trap 
mm and cm particles in the vortex center. Therefore the 
secondary vortices should be able to trap dust particles, 
leading to an asymmetric global dust distribution. 



r (au) 



Figure 10. Pressure perturbation at t = 3500 orbits. The top 
panel shows a cut of the pressure for 0 equals to the secondary 
vortex center. The bottom panel shows a cut of the pressure for 
r equals to the secondary vortex center. The red solid line repre¬ 
sents simulation TRIO, whereas the slate blue dotted-dashed line 
represents simulation IS03MJ. 


6.3. Oph IRS 48 


The system Oph IRS 48 is a good cand idate to host a 
vortex-like structure induced by a planet (|van der Marel 

LmA observa¬ 
tions at U.44 mm revealed a high-contrast asymmetry in 
the disk of this system, which was interpreted as existing 
due to the presen ce of an anticyclonic vortex (van der 
Marel et al. 2013). Besides, this system shows a central 
cavity in CO line observations, which was explained as a 
gap opened by a massi ve planet ( Brown et al.|2Q12 ). 


der Marel et al. (2013) ran a FARGO simulation consid¬ 


ering the parameters of this system to get the gas density 
distribution. Later on the result from the HD simulation 
was used as the initial condition in a dust evolution code 
to get the expected continuum emission. They were able 
to roughly reproduce the ALMA observation; however, 
there is a debate regarding the location of the vortex. If 
the planet is located at 20 AU, the vortex is expected to 
be located at most at ^ 45 AU, nonetheless it is located 
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at ^ 63 AU. 


We run a simula tion using the same setup as van der 


Marel et al. (2013), with the difference that here no vis¬ 
cosity was included, but instead we used thermal relax¬ 
ation, and initialized the disk with i7/r constant. We 
observed vortex formation at the outer edge of the plan¬ 
etary gap, the vortex position is roughly 40 AU; however, 
we did not observe a second generation of vortices. The 
disk considered was much larger (from 2 to 150 AU) than 
the one of our benchmark cases (from 0.25 to 4.0 AU), 
therefore the resolution may not have been sufficient to 
solve the secondary vortex. Figure pT| shows the potential 
vorticity for this simulation after 700 planetary orbits. 
We can see that ( is negative around 60 AU, therefore a 
secondary vortex could be formed in that region. 



r (au) 

Figure 11. Potential vorticity with the Keplerian profile sub¬ 
tracted for t = 700 orbits. The color bar was truncated from 

—0.5 to 0.5 in order to obtain a higher contrast. 

The ratio between the positions of the secondary and 
primary vortices in our benchmark cases is about 1.5. If 
this ratio is fixed, the second generation of vortices in 
the Oph IRS 48 system would be located at ^ 60 AU. 
In order to check whether the secondary vortex was not 
observed due to a numerical resolution problem, we run a 
second simulation considering the same planet-disk setup 
as before, but integrating over a smaller disk size. We 
fixed the inner and outer disk radius, in order to main¬ 
tain the same ratio between the planet orbital distance 
and the boundaries as for the benchmark cases. The 
new simulation has a disk ranging from 5 AU to 80 AU. 
Figure shows the potential vorticity for the smaller 
disk size after 700 planetary orbits. Here, we can see the 
formation of a secondary vortex, located at ^ 62 AU, 
indicating that in the previous case the numerical res¬ 
olution was indeed not sufficient. This new result is a 
promising explanation for the location of the Oph IRS 
48’s vortex, assuming a single planet at ^ 20 AU dis¬ 
tance from the star. It is also important to mention that 
at 700 planetary orbits, the primary and secondary vor¬ 
tices are present. The benchmark cases show that the 
primary vortex get damped before the secondary one, 
thus for later times, just the secondary vortex may be 


present. 



Figure 12. Potential vorticity with the Keplerian profile sub¬ 
tracted for t = 700 orbits. The color bar was truncated from 

—0.5 to 0.5 in order to obtain a higher contrast. 


7. VORTEX LIFETIMES AND BIRTH TIMES 

In this section we obtain the primary vortex lifetime 
and secondary vortex birth time as a function of the ther¬ 
mal relaxation timescale. For this purpose, we first had 
to define when the vortex is born. Once an overdensity 
with negative potential vorticity arises, we declare that 
this overdensity is a vortex. We define as overdensity a 
region that possesses an average surface density at least 
20% higher than the average background surface density. 

The center of the vortex was defined as the position of 
the maximum surface density. We determined the vortex 
edge by looking to the position where the potential vor¬ 
ticity drops to 50% below its value at the vortex center. 
This procedure was done in both radial and azimuthal 
directions. Knowing the dimensions of the vortex, we 
could calculate the average surface density and potential 
vorticity inside the vortex. 

The procedure of finding the vortex and defining its 
border was done until the time that the vortex can not 
be defined as an over density anymore. Once this criterion 
was reached we defined the vortex as dead. In this way, 
the vortex lifetime was defined as the difference between 
the time it is born and the time it dies. 

Figure presents the primary vortex lifetime as a 
function of the different thermal relaxation timescales. 
We obs erve a nonmon otoni c behavior that w as already 
seen b y ~ 


Fu et al. 


(2014) and Les & Lin ( 2QI5[ ). |Les fc Lin 


(2015) explained that the nonmonotonic behavior is due 


to the fact that the vortex lifetime depends on (i) the de¬ 
cay timescale of the RWI, which decreases for increasing 
Ur, and (ii) the vortex growth time, which increases for 
values of Ur up to ^ 5.0 and then decreases for larger 
values. The nonmonotonic nature is a result of this dou¬ 
ble dependence. The double peak, featuring at small 
Ur’s and Ur = 5.0, is due to the nonmonotonic behav¬ 
ior of the vortex growth time for di fferent Ur’s. Higher 
disk temperatures favors the RWI (Li et al. 2000 Lin 
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2012a), nevertheless this effect seems to be important 
just tor larger Ur’s. The dependence of the vortex life¬ 
time with the vortex growth time comes to the fact that 
once the vortex amplitude is very large, it begins to in¬ 
duce shocks, thus the vortex looses energy through shock 
dissipation and starts to decay. The inviscid approxima¬ 
tion may have inffuenced the estimation of the vortex 
lifetimes, si nce it is inversely dependent on the viscosity 


magnitude (|de Val-Borro et al.||2QQ7 Ataiee et ah 


2013 


Fu et ah 2014| ). Nonetheless, it is clear that it did not 


inlluence the qualitatively behavior of vortex lifetimes as 
a function of thermal relaxa tion times c ales, s ince our re¬ 


sult s are i n agreement with Fu et ah (2014) and Les & 


Lin (2015). 
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Figure 13. The lifetime of the primary vortex (red dashed line) 
and the birth time of the secondary vortex (slate blue dotted- 
dashed line) as a function of the thermal relaxation timescale. 


We also plot in Figure the time when the secondary 
vortex is born. A nonmonotonic behavior is also observed 
and the curves are shifted by a few hundreds of planetary 
orbits, with exception for Ur = 10.0. Since the primary 
vortex is born in a scale of tens of planetary orbits, it 
is clear that the secondary vortex is always born before 
the death of the primary vortex, again with exception for 
Ur = 10.0. The time taken for the secondary vortex to 
be born is correlated to the time that the primary vortex 
needs to deplete the mass on its orbit. Therefore it de¬ 
pends on the vortex growth time and the accretion rate 
generated by the primary vortex. This explains the in¬ 
verse dependence of the secondary vortex birth time and 
primary vortex growth time as a function of the thermal 
relaxation time scale. 


8. SUMMARY AND CONCLUSIONS 

Vortices can be formed as a product of the RWI and/or 
radial buoyancy (BI/SBI/CO). The RWI can be trig¬ 
gered in the walls of a planetary gap due to a sharp 
surface density gradient. The disk is buoyantly unsta¬ 
ble when the pressure and entropy gradients have the 
same sign, a thermal relaxation of the order of Ur ^1.0 
also favors vortex amplification. We carried out global 
2D-HD simulations of planet-disk interactions, using the 
PLUTO code. The aim was to study the long-term evo¬ 
lution of planet-induced vortices in inviscid disks and 
initially buoyantly unstable, considering several thermal 
relaxation timescales. Thermal relaxatio n is an impor- 
tant ingredient to sustain radial buoyancy (Petersen et al. 


2QQ7a|b ). It has also a strong impact on amplifying and 


damping vortices. 

We found that radial buoyancy smoothen the surface 
density gradients in the wall of a planetary gap, which 
generates weaker vortices. In this particular physical sce¬ 
nario, radial buoyancy operates against vortex amplifi¬ 
cation and survival. This effect is less pronounced for 
the isothermal and quasi-isothermal states (Ur << 1), 
which is expected, since thermal relaxation is a required 
ingredient to sustain radial buoyancy. The qualitative 
system evolution is similar for different thermal relax¬ 
ation timescales and different planet masses. The major 
difference is regarding the timescales of events (e.g., time 
required for vortex damping and mass transfer). 

The most interesting result from our simulations was 
the formation of a second generation of vortices. The pri¬ 
mary vortex creates an effective a-viscosity that is large 
enough to induce accretion. We obtained a-values in 
the range a = 10“^ — 1Q~^, which agrees with what 

is obtained by the MRI (jDzyurkevich et al.||2010|). The 
accretion process depletes the mass in the primary vor¬ 
tex orbit, creating a density enhancement outwards the 
vortex position. This bump is sufficient to trigger the 
RWI, leading to the secondary vortex formation. This 
result is a promising explanation f or the location of the 


vorte x in the Oph IRS 48 system (van der Marel et al. 
2013), which is located at ^ 63 AU. Previous models 
predicted that the vortex location could be at most at 
^ 45 AU, assuming a single planet at ^ 20 AU. Our 
model suggests that a second generation of vortices can 
be formed at ^ 62 AU, if a massive planet (5Mj) is as¬ 
sumed at 20 AU. We suggest that further works should 
test the formation of a second generation of vortices in 
non-inviscid disks and considering a proper treatment of 
the system’s barycenter location, since these factors may 
influence the generation and sustenance of vortices. 

We observed a nonmonotonic behavior for the vortex 
lifetime as a function of the thermal relaxation times cale. 
This result was al ready observed by Fu et al. (2014) and 


Les & Lin (2015|). The vortex lifetime depends on the 
decay of the RWl and the vortex growth time. The for¬ 
mer decreases as a function of the thermal relaxation 
timescale. The latter increases as a function of the ther¬ 
mal relaxation timescale up to Ur = 5.0, decreasing 
for larger Ur’s. The nonmonotonic behavior and dou¬ 
ble peak observed for the vortex lifetime is a result of 
this double dependence. The birth time of the secondary 
vortex also presents a nonmonotonic behavior. The ap¬ 
pearance of the secondary vortex is correlated to the time 
the primary vortex needs to deplete the mass on its orbit. 
Therefore it is linked to the primary vortex growth time 
and the accretion rate generated by it. It is important to 
remember that we considered an inviscid disk. Previous 
works have shown that the vortex life time is inversely de¬ 


pendent on the viscosity magnitude dde Val-Borro et al 
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m our models is turbulence-triggerea 


All the viscosity 
triggered by the hydrody- 
namical instabilities. The inviscid approximation may 
have quantitatively changed the vortices lifetimes; how¬ 
ever, it did not change the qualitative behavior of vortices 
lifetime as a function of thermal relaxation timescales. 
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